FIELDimageR pipeline:
Image analysis applied to plant breeding

Introduction
- Plant Breeding
- FIELDimageR
FIELDimageR pipeline
- 1) Example 01
- 2) Example 02
Reference

Introduction

Plant Breeding

Images can be used in plant breeding to draw inference about many traits:

Geometric traits (i.e. plant height, leaf area index, lodging, crop canopy cover)
Canopy spectral texture (spectral features)
Physiological traits (i.e., chlorophyll, biomass, pigment content, photosynthesis)
Abiotic/biotic stress indicators (i.e., stomatal conductance, canopy temperature difference, leaf water potential, senescence index)
Nutrients (nitrogen concentration, protein content)
Yield

FIELDimageR

FIELDimageR is a R package to analyze images from plant breeding programs, and allows to:

Crop the image
Remove soil effect
Build vegetation indices
Rotate the image
Build the plot shapefile
Extract information for each plot
Evaluate stand count, canopy percentage, and plant height

FIELDimageR pipeline

1) Example 01

Remote sensing (Potato Breeding)

Jeffrey Endelman and Filipe Matias (UW-Madison)

Orthomosaic using Open Drone Map (ODM)

Follow the OpenDroneMap’s documentation according to your operating system (Windows, macOS or Linux) to install WebODM.

Start WebODM and +Add Project to upload the RGB images from the Tropical Forage Breeding - Embrapa Beef Cattle:

Donwload Example 1: UFV_ODM

After the running process ‘completed’, download the odm_orthophoto.tif and dsm.tif to upload in R. Then follow the pipeline of FIELDimageR.

Data extraction (R/FIELDimageR)

Required packages

FIELDimageR

# Installing 
install.packages("sp")
install.packages("raster")
install.packages("rgdal")
install.packages("readxl")
install.packages("ggplot2")
install.packages("agricolae")
install.packages("reshape2")
install.packages("devtools")
install.packages("lme4")
install.packages("plyr")
devtools::install_github("filipematias23/FIELDimageR")

# Necessary packages 
library(FIELDimageR)
library(raster)
library(readxl)
library(agricolae)
library(reshape2)
library(ggplot2)
library(lme4)
library(plyr)

Donwload: odm_orthophoto.tif


# Uploading file (odm_orthophoto.tif):
EX1<-stack("odm_orthophoto.tif")
plotRGB(EX1)

# Rotating the image with theta 2.3:
EX1.Rotated<-fieldRotate(mosaic = EX1[[-4]],theta = 2.3)

# Removing soil (example using HUE):
EX1.RemSoil<- fieldMask(mosaic = EX1.Rotated)

Evaluating each trial per time:

# Trial: 01
EX1.Crop.1 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic)

Dataset - Field notes - phenotype (Download: ‘EX1_Data.xlsx’)

# Dataset - Field notes - phenotype (.exel or .txt):
Data1<-read_excel('EX1_Data.xlsx',1)

# Field map ID identification (ID="Plot"):
Map1<-fieldMap(fieldPlot = Data1$Plot,fieldColumn = Data1$Column, fieldRow = Data1$Row,decreasing = T)
Map1

# Building the plot shapefile (ncols = 14 and nrows = 10)
x11()
EX1.Shape.1<-fieldShape(mosaic = EX1.Crop.1, ncols = 14, nrows = 10, fieldData = Data1, ID = "Plot", fieldMap = Map1)

# Trial: 02
EX1.Crop.2 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data2<-read_excel('EX1_Data.xlsx',2)
Map2<-fieldMap(fieldPlot = Data2$Plot,fieldColumn = Data2$Column, fieldRow = Data2$Row,decreasing = T)
EX1.Shape.2<-fieldShape(mosaic = EX1.Crop.2, ncols = 14, nrows = 9, fieldData = Data2, ID = "Plot", fieldMap = Map2)

# Trial: 03
EX1.Crop.3 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data3<-read_excel('EX1_Data.xlsx',3)
Map3<-fieldMap(fieldPlot = Data3$Plot,fieldColumn = Data3$Column, fieldRow = Data3$Row,decreasing = T)
EX1.Shape.3<-fieldShape(mosaic = EX1.Crop.3, ncols = 14, nrows = 13, fieldData = Data3, ID = "Plot", fieldMap = Map3)

# Trial: 04
EX1.Crop.4 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data4<-read_excel('EX1_Data.xlsx',4)
Map4<-fieldMap(fieldPlot = Data4$Plot,fieldColumn = Data4$Column, fieldRow = Data4$Row,decreasing = T)
EX1.Shape.4<-fieldShape(mosaic = EX1.Crop.4, ncols = 14, nrows = 10, fieldData = Data4, ID = "Plot", fieldMap = Map4)

# Combining shapefiles
EX1.Shape<-rbind(EX1.Shape.1$fieldShape,EX1.Shape.2$fieldShape,EX1.Shape.3$fieldShape,EX1.Shape.4$fieldShape)
plot(EX1.Shape)

# Building indices ("NGRDI","BGI","GLI","SCI")
EX1.Indices<- fieldIndex(mosaic = EX1.RemSoil$newMosaic,
                         index = c("NGRDI","BGI"), 
                         myIndex = c("(Red-Blue)/Green"))

plot(EX1.Indices$myIndex)
plot(EX1.Shape,add=T)

# Extracting Data
EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI","myIndex")]],
                     fieldShape = EX1.Shape,
                     buffer = -0.05,
                     n.core = 3)
NewData<-EX1.Info$fieldShape@data
NewData

# Making map plots
fieldPlot(fieldShape=EX1.Info$fieldShape,
          fieldAttribute="NGRDI", 
          mosaic=EX1.Indices, 
          color=c("red","blue"), 
          alpha = 0.4,
          round = 2)

Estimation of plant height using a high throughput phenotyping platform. (Donwload: dsm.tif)

# Uploading files from vegetative growth (dsm_1.tif):
DSM <- stack("dsm.tif")
plot(DSM)

# Rotating the image using the same theta (2.3):
DSM.R<-fieldRotate(DSM, theta = 2.3)

# Removing the soil using mask from step 4:
DSM.S <- fieldMask(DSM.R, mask = EX1.RemSoil$mask)

# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(DSM.S$newMosaic, fieldShape = EX1.Info$fieldShape, fun = "mean",n.core = 3)
EPH$plotValue

# Correlation
NewData2<-EPH$fieldShape@data[,c("layer","NGRDI","BGI","myIndex")]
cor1<-correlation(NewData2)
cor1$correlation
cor1$pvalue

# Regression
NewData3<-EPH$fieldShape@data[,c("Trial","Name","Block","Row","Column","layer","NGRDI")]
NewData3$Trial<-as.factor(NewData3$Trial)
NewData3$Name<-as.factor(NewData3$Name)
NewData3$Block<-as.factor(NewData3$Block)
NewData3$Row<-as.factor(NewData3$Row)
NewData3$Column<-as.factor(NewData3$Column)
NewData3$layer<-as.numeric(as.character(NewData3$layer))
NewData3$NGRDI<-as.numeric(as.character(NewData3$NGRDI))

ggplot(NewData3,aes(x=layer,y=NGRDI, col=Trial))+
  facet_wrap(~Trial, scales = "fixed",ncol = 1)+
  geom_point() +
  geom_smooth(method=lm)+
  theme_bw()

fieldPlot(fieldShape=EPH$fieldShape,
          fieldAttribute="layer", 
          mosaic=EX1.Indices, 
          color=c("red","blue"), 
          alpha = 0.4,
          round = 2)

Heritability

# NGRDI

str(NewData3)

mod.T1<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)

mod.T2<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)

mod.T3<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)

mod.T4<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)

NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
                     H2=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))

ggplot(NewData4,aes(x=Trail,y=H2,fill=Trail))+
  geom_bar(stat="identity")+
  geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
  theme_bw()

# Estimated plant height 

mod.T1<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)

mod.T2<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)

mod.T3<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)

mod.T4<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)

NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
                     EPH=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))

ggplot(NewData4,aes(x=Trail,y=EPH,fill=Trail))+
  geom_bar(stat="identity")+
  geom_text(aes(label=round(EPH,2)), vjust=1.6, color="white", size=6)+
  theme_bw()

2) Example 02

Evaluating disease damage (Tomato Breeding)

Fernando Piotto and Jéssica Nogueira (ESALQ/USP)

Donwload: EX2.zip
Pictures of 10 genotypes with Xanthomonas perforans

# Images names (folder directory: "./images/")
pics<-list.files("./EX2/")
pics

# Vegetation indices
index<- c("BI","GLI","NGRDI")

#Choose one image to prepare the pipeline
EX.L1<-stack(paste("./EX2/",pics[1],sep = ""))
EX.L1<-aggregate(EX.L1, fact=4)
plotRGB(EX.L1)

# Shapefile with extent=T (The whole image area will be the shapefile)
EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T)

# Select one index to identify leaves and remove the background
EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T)

plotRGB(EX.L2$newMosaic)

# Select one index to identify demaged area in the leaves  
EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T)

plotRGB(EX.L3$newMosaic)

# Indices
EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index)

# Making a new stack raster with new layers (demage area and indices)
EX.L5<-stack(EX.L4[[index]],(1-EX.L2$mask),(1-EX.L3$mask)) 
names(EX.L5)<-c(index,"Area","Demage")
plot(EX.L5)

# projection=F (Ignore projection. Normally used only with remote sensing images)
EX.L.Info<- fieldInfo(mosaic=EX.L5, 
                      fieldShape=EX.L.Shape$fieldShape, 
                      projection=F)

# Combine information from all images in one table
EX.L.Info$plotValue

* Parallel


# Installing 
install.packages("foreach")
install.packages("parallel")
install.packages("doParallel")

# Necessary packages
library(foreach)
library(parallel)
library(doParallel)

# Number of cores
n.core<-3

# Starting parallel
cl <- makeCluster(n.core, output = "")
registerDoParallel(cl)
system.time({
  EX.Table <- foreach(i = 1:length(pics), .packages = c("raster","FIELDimageR"), 
                               .combine = rbind) %dopar% {
                                 EX.L1<-stack(paste("./EX2/",pics[i],sep = ""))
                                 EX.L1<-aggregate(EX.L1, fact=4)
                                 EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T, plot=F)
                                 EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T, plot=F) 
                                 EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T, plot=F)   
                                 EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index, plot=F)
                                 EX.L5<-stack(EX.L4[[index]],(1-EX.L2$mask),(1-EX.L3$mask)) 
                                 names(EX.L5)<-c(index,"Area","Demage")
                                 EX.L.Info<- fieldInfo(mosaic=EX.L5, fieldShape=EX.L.Shape$fieldShape, projection=F) 
                                 EX.L.Info$plotValue 
                               }})

EX.Table.2<-data.frame(Genotype=do.call(rbind,strsplit(pics,split = ".JPG")),EX.Table[,-1])
EX.Table.2

Dataset: Visual Score
Horsfall-Barratt scale
4 Evaluators
(Download: ‘EX2_Data.txt’)

# Field data
DataEX2<-read.table("EX2_Data.txt",header = T)
DataEX2$Score_HB<-as.numeric(as.character(DataEX2$Score_HB))
DataEX2$Person<-as.factor(DataEX2$Person)
DataEX2$Genotype<-as.factor(DataEX2$Genotype)
DataEX2

ggplot(DataEX2,aes(x=Genotype,y=Score_HB,fill=Person))+
  facet_wrap(~Genotype,scales = "free_x",nrow = 2)+
  geom_bar(stat="identity",position=position_dodge())+
  scale_fill_grey()+
  theme_bw()+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# Regression
DataEX2.1<-ddply(DataEX2,~Genotype,summarise,Score_HB=mean(Score_HB))
DataEX2.2<-merge(DataEX2.1,EX.Table.2,by="Genotype")
DataEX2.2

DataEX2.3<-melt(DataEX2.2,
               value.name = "Index",
               measure.vars = c("Demage","BI","GLI","NGRDI"))
DataEX2.3$Score_HB<-as.numeric(DataEX2.3$Score_HB)
DataEX2.3$Index<-as.numeric(DataEX2.3$Index)
DataEX2.3$variable<-as.factor(DataEX2.3$variable)
DataEX2.3

ggplot(DataEX2.3,aes(x=Index,y=Score_HB,col=variable))+
  facet_wrap(~variable,scales = "free_x")+
  geom_point() +
  geom_smooth(method=lm)+
  theme_bw()

Reference

Matias, F.I. (2019). FIELDimageR Pipeline (tutorial). https://github.com/filipematias23/FIELDimageR

Matias, F.I.; Caraza-Harter, M.; Endelman, J.B. (2020). FIELDimageR: A R Package to Analyze Orthomosaic Images from Agricultural Field Trials. The Plant Phenome Journal. DOI:10.1002/ppj2.20005

FIELDimageR pipeline: Image analysis applied to plant breeding

FIELDimageR pipeline: Image analysis applied to plant breeding

Introduction

Plant Breeding

FIELDimageR

FIELDimageR pipeline

1) Example 01

Remote sensing (Potato Breeding)

Orthomosaic using Open Drone Map (ODM)

Data extraction (R/FIELDimageR)

Required packages

2) Example 02

Evaluating disease damage (Tomato Breeding)

* Parallel

Reference

FIELDimageR pipeline:
Image analysis applied to plant breeding

FIELDimageR pipeline:
Image analysis applied to plant breeding